import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps

A1 = np.loadtxt('200K_25meV.txt')
A2 = np.loadtxt('600K_25meV.txt')
a = np.loadtxt('NWDOS_200K')
b = np.loadtxt('NWDOS_600K')

Area1 = simps(A1[5:68,1],A1[5:68,0])
Area2 = simps(A2[5:65,1],A2[5:65,0])
area1 = simps(a[477:6343,5],a[477:6343,0])
area2 = simps(b[952:12171,5],b[952:12171,0])

print(Area1, Area2,area1,area2)

print(area1/Area1,area2/Area2)

fig, (ax1,ax2) = plt.subplots(2,1,sharex=True,sharey=True)
fig.subplots_adjust(hspace=0)
#left,bottom,width,height = [0.5,0.5,0.35,0.35]
#ax2 = fig.add_axes([left,bottom,width,height])


ax1.errorbar(A1[5:,0],A1[5:,1]/Area1*area1,yerr=A1[5:,2]/Area1*area1,marker='o',markersize=6,label='200K INS',color='blue')
ax1.errorbar(A2[5:,0],A2[5:,1]/Area2*area2,yerr=A2[5:,2]/Area2*area2,marker='s',markersize=6,label='600K INS',color='red')

ax1.fill_between(A1[5:,0],A1[5:,1]/Area1*area1,A2[5:,1]/Area2*area2,color='grey',alpha=0.2,edgecolor="none")
ax1.text(2,0.002,'excess Li DOS',fontsize=20,color='black',rotation=12)

#ax2.plot(a[:,0],a[:,5],label='200K MLMD',color='blue')
#ax2.plot(b[:,0],b[:,5],label='600K MLMD',color='red')
ax2.plot(a[:,0],a[:,5],label=' ',color='blue')
ax2.plot(b[:,0],b[:,5],label=' ',color='red')
ax2.plot(a[:,0],a[:,1],label='200K total, Li DOS',color='blue',zorder=0,ls='dashed')
ax2.plot(b[:,0],b[:,1],label='600K total, Li DOS',color='red',zorder=0,ls='dashed')
ax1.plot(0,b[0,5]*5.0/7.52,label='_nolegend_',color='purple',marker='P',clip_on=False,markersize=15)

print(b[0,5]*5.0/7.52)
ax1.text(0.1,0.0015,'600K QENS',fontsize=15,color='purple')

c = np.interp(a[:,0],b[:,0],b[:,1])
ax2.fill_between(a[:,0],a[:,1],c,color='grey',alpha=0.2,edgecolor="none")


ax1.legend(frameon=False,fontsize=15,loc='upper left',handlelength=1.5)
ax2.legend(frameon=False,fontsize=15,loc='upper left',ncol=2,handlelength=1.5,columnspacing=-0.5)

ax1.set_xlim(0,4)
ax1.set_ylim(0,0.01)
ax1.set_yticks([0.00,0.004,0.008])
ax1.tick_params(which='both',direction='in',labelsize=20)

ax2.set_xticks([0,1,2,3,4])
ax2.tick_params(which='both',direction='in',labelsize=20)
ax2.minorticks_on()

fig.add_subplot(111, frameon=False)
# hide tick and tick label of the big axis
plt.tick_params(labelcolor='none', which='both', top=False, bottom=False, left=False, right=False)
plt.tick_params(axis='y',pad=50)
plt.xlabel('Energy (meV)',fontsize=20)
plt.ylabel(r'DOS (meV$^{-1}$)',fontsize=20)

plt.savefig('NWDOS_zoom_v6.png',format='png',bbox_inches='tight',dpi=600)
#plt.savefig('NWDOS_zoom_v6.pdf',format='pdf',bbox_inches='tight',dpi=600)
plt.savefig('NWDOS_zoom_v6.eps',format='eps',bbox_inches='tight',dpi=600)
plt.show()
